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ABSTRACT 


A  block-structured  solution  scheme  is  developed  for  the  analysis  of 
three-dimensional  transonic  flows.  The  scheme  is  based  on  the  solution 
of  potential  flow  equations  for  individual  blocks  representing  part  of 
the  flow  field.  Based  on  a  previously  developed  block-structured  grid 
generation  scheme,  appropriate  computational  grids  are  generated  for  each 
of  the  blocks  depending  on  the  complexity  of  the  local  flow  field.  The 
equations  are  then  solved  to  provide  a  solution  of  a  large  problem  in 
terms  of  an  assembly  of  smaller  problems  for  each  block. 

Numerical  results  illustrate  the  applicability  of  the  method  for  a 
three-dimensional  flow  field  around  a  wing  profile  (NACA0012).  Different 
block  structures  are  analyzed  to  demonstrate  the  robustness  and  the  ac¬ 
curacy  of  the  developed  method.  Finally  a  three-dimensional  wing-body 
configuration  is  analyzed  and  the  results  are  compared  with  previously 
obtained  single  block  solutions. 

The  method  is  expandable  to  the  solution  of  Euler  and  Navier-Stokes 
equations.  It  is  also  suited  to  be  executed  in  a  parallel  processing 
environment. 
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I.  INTRODUCTION 


The  solution  of  three-dimensional,  transonic  flows  around  complex 
aircraft  configurations  requires  considerable  computational  effort.  The 
capability  to  solve  such  large  problems  relies  heavily  on  the  availability 
of  larger  and  faster  computers.  Over  the  last  twenty  years,  considerable 
progress  has  been  made  in  terms  of  hardware  available  for  solving  compu¬ 
tational  fluid  dynamics  problems.  Rakich  [1],  in  his  introduction  to  the 
1973  AIAA  Computational  Fluid  Dynamics  Conference,  compared  machines  like 
IBM  360/67  to  the  kinds  of  CDC  Star  and  Burroughs  ILLIAC.  At  the  time, 
parallel  and  vector  machines  were  recently  introduced  and  being  applied 
to  the  solution  of  computational  fluid  dynamics  problems.  In  1983, 
further  hardware  development  has  resulted  in  vector  machines  like  CRAY 
X  MP  or  CYBER-205  [2].  Such  comparisons  may  indicate  the  kind  of  progress 
in  computer  power  we  may  expect  in  coming  years.  It  snould  also  be  men¬ 
tioned  that  these  types  of  machines  have  also  become  much  more  widely 
available  to  the  researchers  in  the  field.  In  the  next  five  years,  we 
expect  the  major  emphasis  in  hardware  development  to  be  in  terms  of  multi¬ 
processing  capabilities  and  larger  memories  rather  than  drastic  changes 
in  computational  speeds. 

The  computation  of  flow  problems  using  vector  machines  has  resulted 
in  comparisons  of  computational  speeds  with  other  machines.  For  example, 

Chima  and  Johnson  [3],  compared  Euler  and  Navier-Stokes  solutions  for 
transonic  flows  through  a  cascade  of  airfoils  on  two  different  computers: 

IBM  370/3033  and  CRAY  1-S.  In  their  comparisons,  they  have  employed  different 
levels  of  grid  refinement  for  which  they  employed  a  Lax-Wendroff  scheme. 

It  can  be  seen  from  such  comparisons,  however,  that  the  evaluation  of  the 
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efficiency  of  a  particular  hardware  configuration  is  not  straight¬ 
forward.  It  depends  on  the  computational  scheme  and  it  is  also  strongly 
related  to  the  computational  grid  employed  in  the  analysis.  It  is  pos¬ 
sible  to  define  a  computational  fluid  dynamics  problem  on  a  box-type  grid 
structure  with  MxNxK  grid  points  in  each  direction.  One  can  then  specify 
the  efficiency  of  a  numerical  scheme  for  such  a  rectangular  grid  struc¬ 
ture  and  for  a  given  hardware  configuration. 

A  second  problem  to  be  considered  is  the  number  of  grid  points  in  a 
single  block  which  can  be  fitted  in  a  particular  vector  machine.  If  the 
strategy  is  to  fit  a  single  block  with  a  regular  grid  structure  around  an 
aircraft  geometry,  the  number  of  grid  points  required  to  model  the  complete 
flow  field  can  grow  very  rapidly.  Attempts  are  being  made  to  improve  this 
situation  by  embedding  blocks  with  several  levels  of  grid  refinement  in¬ 
side  a  regular  but  coarser  grid  structure  [4].  Until  now,  the  most 
popular  strategy  for  solving  large  computational  fluid  dynamics  problems 
has  been,  a)  to  develop  computational  schemes  which  are  fast  on  regular 
grids  and,  b)  to  model  irregular  geometries  using  regular  grid  structures. 

Our  work  has  been  directed  in  applying  finite  element  methods  for 
solving  fluid  mechanics  problems.  This  method  provides  an  alternate  ap¬ 
proach  over  the  ones  discussed  above  for  the  solution  of  the  same  problem. 
One  of  the  important  features  of  this  method  is  the  possibility  of  using 
irregular  grids  that  can  be  fitted  efficiently  to  simulate  irregular  flow 
fields  around  complex  geometries.  It  provides  sufficient  flexibility  to 
design  a  computational  grid  around  a  complex  aircraft  geometry  [5].  It 
also  provides  the  freedom  to  modify  an  existing  grid  to  provide  a  better 
approximation  basis  for  the  employed  numerical  scheme. 

The  applications  of  irregular  grids,  however,  require  development  of 
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numerical  schemes  where  no  regularity  of  the  nodal  connectivities  is 
assumed.  Our  work  has  been  aimed  at  developing  such  numerical  schemes 
[6]  where  we  can  exploit  the  advantages  of  irregular  yet  more  efficient 
grids.  As  it  will  be  discussed  in  the  following  sections,  rather  than 
employing  point  or  line  relaxation  schemes,  we  use  block-relaxation 
schemes.  Such  schemes  produce  a  uniform  convergence  rate  for  all  of  the 
eigenvalues  of  a  differential  operator  inside  a  block.  It  is  not  affected 
by  the  occurence  of  an  irregular  grid  structure  inside  a  block.  Such 
schemes  have  two  basic  constraints  for  solving  large  problems:  a)  the 
size  of  the  block  operator  grows  proportional  to  the  size  of  the  block 
and,  b)  the  geometric  definition  of  an  irregular  grid  requires  much  more 
data  than  a  regular  grid,  thus  considerable  information  has  to  be  generated 
or  read-in  during  each  relaxation  step. 

Based  on  the  above  considerations,  we  have  been  working  on  develop¬ 
ing  a  block-by-block  solution  scheme.  The  main  objective  of  this 
approach  is  to  divide  a  large  problem  into  smaller  components  in  terms 
of  a  series  of  blocks.  Rather  than  attempting  to  solve  a  large  problem 
most  efficiently  on  a  single  processor,  we  divide  the  problem  into  smaller 
ones  and  try  to  develop  an  "intelligent"  strategy  which  is  suitable  for 
parallel  processing.  Parallel  processing  for  large  systems  is  a  popular 
subject  addressed  by  many  researchers  today.  Our  main  objective  is  to 
exploit  the  physical  characteristics  of  the  problem  where  each  block 
corresponds  to  a  sub-volume  in  physical  space.  We  developed  a  "sub¬ 
structuring"  scheme  where  each  of  the  "sub-structures"  corresponds  to  a 
particular  flow  region.  One  can  then  design  grids,  use  efficient  solu¬ 
tion  schemes  for  each  of  the  sub-regions  depending  on  the  characteristics 
of  the  flow  field,  and  allocate  computer  resources  in  a  parallel  process- 
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ing  environment  in  a  most  efficient  manner.  We  can  summarize  our  progress 
in  this  area  along  the  following  lines: 

•  Development  of  a  block-structured  grid  generation  scheme  [5], 

•  Development  of  a  variational  formulation  for  solving  potential, 
Euler,  and  Navier-Stokes  equations  in  a  unified  form  [7],  and 

•  Development  of  a  block-structured  solution  scheme  for  solving  the 
steady-state  flow  equations  through  a  relaxation  scheme. 

In  this  report  we  will  present  only  a  short  summary  of  the  first  two  steps. 
We  will  mostly  concentrate  on  the  basic  features  of  the  block-structured 
solution  scheme  and  its  application  to  three-dimensional  transonic  poten¬ 
tial  flows. 

A.  Block-Structured  Grid  Generation  Scheme 

The  generation  of  grids  around  a  complex  aircraft  configuration  is 
not  a  straight-forward  task.  It  requires  considerable  understanding  of 
the  surrounding  flow  field.  One  has  to  provide  a  grid  which  will  produce 
accurate  results  at  all  critical  regions.  Yet,  excessive  refinement  of 
three-dimensional  grids  is  not  possible  even  when  the  largest  computers 
are  available.  One  has  to  be  able  to  generate  a  computational  grid  in 
such  a  manner  that  a  block  of  grid  points  can  be  controlled  and  modified 
after  inspection.  The  optimum  grid  configuration  for  each  of  the  critical 
flow  regions  can  be  quite  different,  yet,  they  have  to  be  connected  to 
each  other.  The  grid  generation  techniques  to  be  employed  for  this  purpose 
should  be  sufficiently  general  to  model  any  complex  configuration. 

In  terms  of  generating  a  block-structured  grid  generation  scheme,  a 
finite  element  approach  provides  certain  advantages.  As  discussed  above, 
since  irregular  grids  are  allowed  in  a  finite  element  solution  scheme,  it 


is  easier  to  design  grids  individually  for  each  block  and  then  assemble 
them  together.  In  comparison,  the  global  mapping  techniques  used  for  grid 
generation  requires  additional  computations  to  provide  appropriate  coupling 
of  the  blocks  [8]. 

The  developed  block-structured  grid  generation  scheme  can  be  summariz¬ 
ed  in  terms  of  the  following  steps:  If  the  aircraft  can  be  considered  as 
as  assembly  of  the  body,  engines,  wings,  and  many  other  components,  one 
has  to  provide  appropriate  grids  for  each  of  these  components.  Figure  1 
illustrates  the  work  performed  for  generating  the  details  of  the  air¬ 
craft  geometry.  As  can  be  seen  in  this  figure,  several  sections  are 
digitized  to  obtain  the  geometry  and  a  three-dimensional  representation  of 
the  body  is  obtained.  A  series  of  blocks  are  then  constructed  around  the 
aircraft  to  model  the  flow  field.  The  block  structure  for  this  problem 
is  shown  in  Figure  2.  As  can  be  seen  from  this  figure,  the  block  struc¬ 
ture  is  irregular,  i.e.,  it  includes  openings  between  the  blocks,  it  has 
irregular  blocks  (tetrahedrons)  and  finally  some  of  the  blocks  are  voids. 
The  developed  grid  generation  scheme  can  treat  such  an  irregular  block 
structure.  One  can  design  quite  irregular  grids  in  each  of  the  blocks, 
yet  attach  the  blocks  in  an  organized  fashion,  as  shown  in  Figure  3. 

Details  of  this  approach  can  be  found  in  reference  [5]. 

The  main  objective  of  the  work  described  in  the  present  report  is  to 
determine  the  optimum  computational  procedure  in  terms  of  number  of  compu¬ 
tations  and  number  of  data  transfers  for  analyzing  such  irregular  block 
structures.  Generation  of  a  block-structured  grid  provides  the  most 
natural  way  to  develop  a  block-structured  solution  scheme.  It  also  pro¬ 
vides  an  insight  for  designing  the  most  efficient  solution  scheme  for  a 
particular  problem. 
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QUADRATIC 


A  typical  Block 


)  One  surface  collapsed.  b)  Two  surfaces  collapsed 


Figure  3.  Capabilities  of  grid  generation  scheme  for 
modeling  irregular  block  structures. 


d)  Coupled-surface  feature  for 
treating  radial-type  grids. 


Figure  3.  'Continued. 


VOID  BLOCK 


B.  Solution  of  Flow  Problems  in  Terms  of  Clebsch  Variables 

The  ultimate  objective  in  computational  fluid  mechanics  is  to  solve 
Navier-Stokes  equations  for  three-dimensional  turbulent  flows.  Since 
such  a  proposal  is  computationally  costly,  a  general  and  yet  quite  effi¬ 
cient  approach  is  to  solve  these  equations  only  in  flow  regions  where 
viscous  effects  are  important.  One  of  the  most  successful  applications 
of  computational  fluid  mechanics  has  been  in  the  area  of  solving  potential 
flow  problems  together  with  boundary  layer  corrections.  We  have  attempted 
to  achieve  such  a  simplification  in  a  general  manner  by  using  a  new  formu¬ 
lation  of  the  physical  problem.  Rather  than  using  the  velocities  and  pres¬ 
sure  as  the  primitive  variables,  we  decomposed  the  velocity  vector  into  ir- 
rotational,  rotational-inviscid  and  viscous  components  by  using  a  Clebsch 
transformation  [9].  This  enables  the  solution  of  the  general  problem  in 
terms  of  potential,  Euler,  and  Navier-Stokes  approximations  at  different 
flow  regions  which  in  our  case  are  defined  as  blocks.  The  objective  is  to 
provide  general  yet  efficient  solution  schemes  as  well  as  appropriate  com¬ 
putational  grids  in  a  block-structured  flow  domain  for  all  three  cases. 

For  analyzing  potential  flows,  it  is  common  to  employ  the  velocity 
potential  as  the  primitive  variable.  One  can  write  the  velocity  vector 
in  terms  of  a  velocity  potential  as  follows: 

u  =  v<j>  (1) 

and  proceed  to  solve  the  conservation  of  mass  equation.  In  this  case, 
the  conservation  of  momentum  and  energy  are  automatically  satisfied  [7]. 

A  second  step  in  the  approximate  solution  of  flow  problems  is  the  Euler 
equations.  In  this  case,  one  can  write  the  velocity  field  in  terms  of 
the  entropy  and  total  enthalpy  as  follows: 

u  =  +  Hvn  +  Svx  (2) 
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where  <j>,  n,  A  are  a  set  of  Lagrange  multipliers  [9].  Here  the  equations 
for  the  conservation  of  mass,  entropy  and  total  enthalpy  have  to  be  solved. 
Finally,  for  viscous  flows,  the  velocity  field  can  be  written  in  terms  of: 

u  =  +  HVn  +  Svx  +  viscous  terms.  (3) 

Here  the  viscous  dissipation  has  to  be  calculated  together  with  the 
shear  stresses  on  the  solid  walls  for  a  given  velocity  field. 

In  the  foregoing  three  levels  of  approximation,  one  can  observe 
that  the  number  of  variables  and  the  number  of  equations  increases 
parallel  to  the  complexity  of  the  model.  If  we  assume  that  in  one  flow 
region  (a  block),  we  will  be  solving  Euler  equations, while  at  a  neighbor¬ 
ing  block  the  potential  flow  equations  are  sufficiently  accurate,  the  inter 
face  boundary  requires  specification  of  S  and  H  to  be  a  constant  along 
the  boundary.  Only  one  equation  is  solved  in  the  irrotational  flow  region, 
while  three  equations  are  solved  in  the  block  with  the  rotational  flow. 

Another  point  to  be  remembered  is  the  relationship  between  the  flow 
models  and  the  computational  grids.  One  can  easily  realize  that  more 
complicated  models  require  much  finer  grids.  For  example,  the  solution 
of  two-dimensional,  transonic  potential  flows  around  a  semi-circle  requires 
a  much  coarser  grid  than  the  one  for  solving  the  Euler  equations  for  the 
same  geometry  [9].  If  the  convection  of  vorticities  generated  by  the  shock 
is  considered,  one  has  to  design  much  finer  grids  to  capture  such  a 
phenomenon. 

C.  The  Necessity  for  Block-Structured  Solution  Schemes 

For  the  development  of  a  block-structured  scheme  the  methodology 
described  above  was  employed.  The  computational  strategy  is  based  on  the 
capability  of  generating  grids  with  varying  degrees  of  refinement  for  each 
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block  and  solve  potential,  Euler,  or  Navier-Stokes  equations  for  each 
block  depending  on  the  characteristics  of  the  flow  field.  As  it  will  be 
discussed  later,  accurate  and  efficient  implementation  of  the  boundary 
conditions  between  neighboring  blocks  is  extremely  important  in  a  block- 
structured  scheme. 

The  most  important  consideration  which  necessitates  the  use  of  a 
block-structured  solution  scheme  is  the  size  of  the  problem.  We  do  not 
expect  to  have  computers  large  enough  to  fit  the  problems  we  are  interested 
in  into  their  main  memory  during  the  next  five  years.  In  terms  of  storing 
the  geometry  information  describing  an  irregular  geometry  in  a  general 
form,  data  transfer  becomes  a  major  problem.  If  one  employs  auxiliary 
storage  to  store  large  amounts  of  data,  the  computational  scheme  can  easily 
become  I/O  bound  when  large  amounts  of  data  have  to  be  transferred. 

Depending  on  the  size  of  the  problem  and  the  type  of  information  to  be 
stored,  one  can  even  sometimes  deplete  these  types  of  storage  facilities. 

The  approach  we  have  taken  here  is  to  develop  a  block-structured 
solution  scheme,  where  it  is  possible  to  use  the  available  computer 
resources  in  a  most  efficient  manner.  For  a  given  problem  one  has  to 
make  decisions  in  terms  of  designing  a  computational  grid.  Depending  on 
the  size  of  the  problem,  the  computer  resources  available,  and  the  char¬ 
acteristics  of  the  flow  field,  we  can  decide  on  the  size  of  the  blocks. 

In  the  developed  numerical  procedure,  based  on  physical  insight,  one  can 
decide  the  manner  in  which  iterations  should  be  performed  and  the  data  should 
be  transferred  between  the  blocks.  If  we  have  a  supersonic  pocket  develop¬ 
ing  inside  a  single  block,  we  may  choose  a  particular  iterative  procedure  for 
solving  this  problem.  We  may  decide  to  iterate  on  that  block  more  often  than 
the  others.  The  objective  here  is  to  develop  "intelligent"  schemes  where 


efficiency  can  be  improved  as  we  learn  more  about  the  details  of  the 
fluid  mechanics  problem  and  we  can  plan  and  revise  our  computational 
strategy. 

In  this  report  we  will  discuss  only  the  solution  of  the  potential 
equation  using  a  block-structured  solution  scheme.  However,  we  will 
comment  on  its  generality  and  its  applications  to  more  complex  flow 
models. 


II.  A  BLOCK-STRUCTURED  SOLUTION  SCHEME  FOR  POTENTIAL  FLOWS 


A.  Potential  Flow  Problem 

The  analysis  of  compressible,  irrotational  flows  requires  only  the 
solution  of  the  following  conservation  of  mass  equation: 

v(pu)  =0  in  ft  ,  (4) 

where  u  is  the  velocity  vector,  p  is  the  density  and  ft  is  the  flow  region 
to  be  analyzed.  The  density  is  a  function  of  the  local  velocity  which  can 
be  written  as  follows: 

1 

p  =  C(K2  -  u.u)Y_1  (5) 

In  the  above  equation,  C,  K  and  y  are  known  constants.  The  boundary  condi¬ 
tions  require  the  specification  of  the  normal  mass  flux  on  the  boundary 
surface  r. 

pu*n  =  f  on  r  (6) 

where  f  is  a  known  function  specified  on  the  boundary  and  n  is  the  unit 
normal  vector.  By  using  the  condition  for  irrotationality,  one  can  sub¬ 
stitute  the  velocity  potential  to  eliminate  the  velocity  vector  as  follows: 

u  =  vj  (7) 

The  conservation  of  mass  equation  then  becomes  second-order, 

v-(pv$)  =0  in  ft  (8) 

where  the  velocity  potential  has  to  be  assigned  an  arbitrary  value  at  one 
point  to  remove  the  singularity  due  to  the  introduction  of  equation  (7). 
Since,  the  conservation  of  mass  equation  has  to  be  satisfied  for  the  entire 
flow  domain,  over  a  closed  boundary,  i.e., 

<j>r  pu*n  dr  =  0  (9) 

the  specified  value  of  <t>  can  be  arbitrary. 
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The  objective  of  the  present  analysis  is  to  develop  a  relaxation 
scheme  for  obtaining  a  solution  to  the  above  steady-flow  equations. 

B.  Formulation  of  the  Problem  in  a  Block-Structured  Form 

For  developing  a  block-structured  solution  scheme,  we  will  define  the 
problem  over  a  series  of  flow  regions  which  we  will  call  blocks  in  the 
following  form: 

V«  ( pv$.j )  =  0  in  ft.. 

pn-v$i  =  f.  on  r.  (10) 

when  n.  and  r\  are  the  block  volume  and  block  surface  for  each  block, 
where  i  indicates  the  block  number.  The  boundaries  can  be  further  classi¬ 
fied  in  two  parts:  the  boundaries  between  the  neighboring  blocks  and  the 
global  boundary  of  the  entire  flow  region.  We  can  distinguish  the  corre¬ 
sponding  boundary  conditions  as  follows: 


pn*v$>..  =  f . 

on  r..0 

(11) 

pn*v$i  =  g.. 

on  r.. c 

(12) 

where  r..0  is  the  global  boundary  and  r^C  is  the  inter-block  boundary  for 
the  i^*1  block.  Of  course,  some  blocks  may  have  no  connection  to  global 
boundaries. 

The  problem  can  then  be  defined  as  the  solution  of  conservation  of 
mass  equation  (10)  for  each  block,  together  with  the  determination  of  un¬ 
known  inter-block  boundary  fluxes  g. .  This  can  be  achieved  through  an 
additional  constraint  which  specifies  that  the  velocity  potentials  are 
continuous  across  the  inter-block  boundaries.  We  can  define  the  corre¬ 
sponding  variational  problem  in  the  following  form: 

fiif  =  J  /v  pv^-vs^.  dV  +  JT o  f.<si>.  dr  +  /p c  g i 5<t-i  dr  =  0  (13 
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By  taking  variations  with  respect  to  $,  we  can  derive  the  differential 
equation  in  equation  (9)  and  the  boundary  conditions  in  equations  (11) 
and  (12).  Also,  since  we  assumed  that  g.  is  an  unknown,  we  can  take  a 
variation  with  respect  to  g^ .  This  produces  the  additional  constraint 
equation  in  terms  of  the  compatibility  of  velocity  potentials  across 
neighboring  block  boundaries  as  follows: 

I  /rc  sg.  dr  =  0  (14) 

Since,  the  same  g^  is  employed  for  neighboring  blocks,  this  constraint 
will  involve  the  velocity  potentials  which  correspond  to  the  grid  points 
on  such  blade  pairs  as  follows: 

-1,j  "  -n,j  =  0  (,5) 

where  i  and  n  are  the  two  neighboring  blocks  and  j  is  the  number  for  the 
boundary  surface  between  these  two  blocks. 

At  this  point  one  can  introduce  the  finite  element  formulation,  where 
the  flow  region  is  divided  into  first  a  series  of  blocks  and  then  a  group 
of  finite  elements.  The  distribution  of  velocity  potential  over  each 
element  for  block  number  (i)  is  approximated  by: 

$i(x,y,z)  =  Nk(x,y,z)  *1k  (16) 

where  Nk  is  called  a  shape  function  (a  simple  polynomial)  and  ^  is  the 
nodal  value  of  the  velocity  potential  at  node  number  k.  One  can  also  ap¬ 
proximate  the  unknown  boundary  fluxes  along  the  block  interfaces  in  terms 
of  nodal  flux  values  as  follows: 

9j(s,t)  =  Nk(s,t)  Xjk  (17) 

where  Nk  is  a  two-dimensional  shape  function  for  node  k  and  (s,t)  are  the 
two-dimensional  coordinates  along  the  boundary  surfaces.  After  substituting 


the  finite  element  approximations  and  taking  variation  with  respect  to 
nodal  values  of  velocity  potentials,  one  can  write  the  following  set  of 
discrete  equations 

Ji  u  +  sf,j  ij  *  -Fi  <18> 

Where  A,  F,  C»  $  *  are  global  matrices  for  each  block  (i)  defined  by 

assembling  element  matrices  for  each  element  e. 


-i,e  =  Al.  p(-'x-’x  +  -’y^’y  +  dV 


Ei,e=/r?  pefe  -  dr  <20> 

C..  i  e  -  -  /rc  pJ  N*  dr  (21) 

i,J,e  1  i  ,e  e 

-  [Xk]j  (22) 

Here,  i  indicates  the  number  of  the  block,  j  indicates  the  number  of  the 
sui  face,  k  is  the  node  number  and  e  is  the  element  number.  For  p>o,  the 
coefficient  matrix  in  equation  (19)  is  symmetric  and  positive-definite. 

The  variation  with  respect  to  nodal  values  of  boundary  fluxes  Uk), 
produces  a  set  of  constraint  equations  defining  the  compatibility  of  velo¬ 
city  potentials  between  the  blocks: 

(23) 

As  can  be  seen  from  the  above  equation,  matrix  C.  .  picks  up  all  the  nodes 

■  *J 

on  surface  j  located  on  the  block  i  and  assembles  the  constraint  conditions. 
In  coupled  form  equations  (17)  and  (20)  can  be  written  as  follows: 


W  eij]  [♦,]  Ci) 


,rcu]  to]  c»j] 
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For  a  block  structure  shown  in  Figure  4,  the  equation  (24)  can  be  written 
in  the  following  form: 


A,  0  g  o 

i  — 

.o, 

Q  a2  g  c2j  z 

-2 

-2 

?  -  -3  g  -3,2 

-*3 

= 

h 

£l,l  -2,1  -  -  - 

*1 

o 

-  -2,2  -3,2  -  - 

-2 

o 

(25) 


Here,  there  are  three-blocks  and  two  inter-connecting  surfaces.  In  the 

above  equation,  the  coupling  matrices  for  each  row  of  A  equations  represent 

the  contribution  to  a  surface  from  the  nodes  located  on  the  two  neighboring 

blocks.  As  indicated  above,  unknowns  are  the  nodal  velocity  potentials 

<j>k  and  the  nodal  boundary  fluxes  A^.  (^  i  =  1,3)  are  the  potential 

vectors  for  each  block  and  (A-  j  =  1 ,2)  are  the  vectors  for  each  surface. 

“J  * 

In  general,  the  surfaces  can  be  described  by  an  independent  set  of 
nodes  and  the  nodal  values  of  A.  vectors  do  not  have  to  be  necessarily 

J 

positioned  on  the  same  physical  point  in  space  as  of  the  elements  of  the 
$  vectors.  However,  in  the  applications  presented  in  this  paper,  the  same 
grid  point  locations  were  used  for  both  $  and  A  vectors. 

C.  Basic  Considerations  for  the  Development  of  an  Iterative  Solution 
Scheme  for  Block-Structured  Equations 

One  can,  of  course,  attempt  to  solve  the  coupled  equations  directly 
in  the  form  of  equation  (24).  Matrices  A.,  correspond  to  individual  co¬ 
efficient  matrices  for  each  block.  Most  of  the  computational  effort  in 
solving  such  a  system  will  involve  the  treatment  of  the  coupling  matrices 
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between  the  blocks.  In  the  developed  procedure,  we  attempt  to  solve  the 
coupled  equations  by  using  an  iterative  procedure.  For  the  description 
of  this  iterative  procedure,  let  us  consider  a  simpler  system  consisting 
of  two  blocks  and  a  single  coupling  surface.  In  this  case,  the  equations 
can  be  written  as  follows: 


-l  9  9i,i 

*i 

"-i" 

Q  a2 

-2 

II 

-2 

-1,1  -2,1  - 

1  1 

1  * 

'  io 

1 

(26) 


Here  A-j ,  and  represent  the  Laplace  operator  for  the  two  blocks.  C-j  ^ 
and  i  are  matrices  for  defining  the  coupling  between  the  two  block'' 

Let  us  also  consider  a  linearized  discrete  variational  problem  in  the 
following  form: 

n  =  \  $lMl  +  \  -2-2-2  +  ^-1,1  $ 1  +  -2,1  -2* 

-  -  *2-2  (27) 

It  can  be  shown  that  the  variation  of  the  above  functional  with  respect  to 

<j>1 ,  $2  and  x  produces  the  above  set  of  equations. 

At  this  point  we  can  define  a  second  variational  problem  as  follows: 

★ 

Instead  of  the  unknown  flux  vector  x,  we  assume  a  known  flux  vector  X  . 
Then  we  can  write  the  following  variational  problem: 

n  =  1  Ml  +  \  -2  -2-2  +  -1  ^-1 ,1-2  +  -2,1-2^ 

>  I  *t 
"  ’  I  : 1  "  *2  -2 

*  * 

Taking  variation  wi  i  respect  to  <t>,  and  <J>9,  we  can  write: 


(28) 


(29) 


-1*1  +  ^U1A1  =  -1 


-2-2  +  -2,1-1  =  -2 


We  can  prove  a  certain  Theorem,  from  the  solution  of  the  above  two  varia¬ 
tional  problems: 

★ 

Theorem.  If  X  is  any  arbitrary  flux  distribution  on  a  boundary 
surface,  the  solution  of  equation  (29)  produces  upper  and  lower  bounds  to 
the  solution  of  the  original  problem  of  equation  (26). 

The  validity  of  the  above  statement  can  now  be  shown.  The  solution 
of  the  above  equations  can  be  written  as: 


*  -1  -It* 

*1  =  -1  -1  “  -1  -1,1*1 


*  -1  -it* 

-2  =  -2  -2  "  -2  -2,1*1 


(30) 


The  nodal  potential  vectors  which  correspond  to  the  nodes  on  the  interface 
can  be  written  from  the  above  solutions  as  follows: 


Si  ,1*1  =  -1 ,1-1  ^-1  '  -1 ,1  -1 1  -1 ,1  - 
-2,1-2  =  -2.1-21  -2  "  ^2, 1-2^2, 1* 

We  can  write  the  exact  solution  of  equation  (26)  in  a  similar  form, 

Si  ,1*1  =  Si, 1^1  - 

-2 ,1-1  =  -2,1 -21  ^2  '  h, 1^2, 1*1 


(31) 


(32) 


The  difference  of  equations  (31)  and  (32)  produce  the  following  rela¬ 
tionship: 

l,t 


Si  ,1^1  "  *1^  "  ?i  j^i  9]  '  A  ) 


lrt 


-2  1^-2  ~  -2^  =  “2  1-2  -2  1  ^ -  -  ) 


(33) 
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If  we  assume  that  ^  and  $2  are  close  to  the  exact  solution  and  can  be 

written  in  the  following  form: 

* 

-1  =  h-1  +  sma^  terms 
★ 

-2  =  ^2-2  +  sma^  terms  (2 

where  and  k2  are  unknown  constants,  equation  (33)  can  be  re-written 
as  follows: 


(kl  '  ^  -1,1  -1  =  -1,1-1  ^1,1  ~  h  ) 

( k2  ~  1)  C2  i  $2  -  C2  -jA2  C2  i  ( ^  -  A  )  (35) 

By  using  the  compatability  condition  for  the  exact  solution  from  equation 

(26)  as  follows: 

9]  i  ^1  +  ^2  1  -2  =  0  *  (36) 

and  since  both  matrices  A-|  and  A2  are  positive  definite  and  by  using  equa¬ 
tion  (35)  we  can  show  that  when 


k-j  >  1 


k2  <  1 


k1  <  1 


k2  >  1 


i.e.,  $.|  and  $2  provide  upper  and  lower  bounds  to  the  solution  of  velocity 
potentials  on  the  boundary. 

A  similar  proof  can  be  made  for  a  second  problem.  This  time,  let  us 

assume  that  we  start  with  a  set  of  approximate  velocity  potential  vectors 
★  ★ 

for  each  block,  and  $2>  such  that  the  compatability  condition  at  the 
boundaries  is  satisfied  as  follows: 

C-j  i  $.j  +  C2  i  $2  =  0  (38) 
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which  means  that  for  the  approximate  solution,  the  nodes  at  both  sides 
of  the  interface  between  two  blocks  possess  the  same  $  distribution, 
yet  equation  (32)  is  not  satisfied.  In  this  case,  by  solving  equation 
(38)  together  with  the  following  set  of  equations. 


-1-1  +  ^1, 1-1,1  =  -1 
-2-2  +  -2, 1-1, 2  =  -2 


(39) 


we  obtain  a  different  set  of  flux  vectors  X-j  -|  and  ^  on  each  side 
the  boundary  surface.  If  we  solve  equation  (39)  for  the  velocity  poten¬ 
tial  vector  as 


*1  =  -/-I  '  -1^1, 1-1,1 

*  -1  -It* 

-2  =  -2  -2  "  -2  -2,1-1 ,2 


(40) 


and  substitute  into  equation  (38),  we  can  write  the  following  relationship: 

Sl.l^Ei  "  -1 ,1-1 1  -1 ,1-1 ,1  +  -2,1 -2^ -2  "  -2,1-2^ -2,1-1 ,2  =  0  (41) 

On  the  other  hand,  the  summation  of  the  two  equations  in  equation  (32) 
produces  the  following  relationship 

-1,1-1  +  -2,1-2  =  -lJ-Z-l  '  -1 ,1  -1 1  -l" ,1-1  +  -2^-2  -2 


’  -2,1-21-2,1-1 ,1 

Combining  equations  (41)  and  (42)  one  can  then  write, 

-1 , 1  -1 1  ,  1  ^1  "  ^1,1*  =  "  -2,l-2]-2,l  ^-1  '  -1,2^ 


(42) 


(43) 


Again  by  making  a  similar  approximation  as  in  the  case  of  equation  (35) 

and  by  recognizing  that  both  matrices  and  ^  are  positive-definite, 

*  * 

one  can  show  that  vectors  x-j  -j  and  X-j  ^  produce  upper  and  lower  bounds  to 
the  correct  flux  vector  at  the  boundaries  which  can  be  stated  as  follows: 
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★ 

Theorem.  If  $  is  any  arbitrary  velocity  potential  distribution  on 
a  boundary  surface,  the  solution  of  equation  (39)  produces  upper  and  lower 
bounds  to  the  solution  of  the  original  problem  of  equation  (26). 

8ased  on  the  above  discussions,  an  iteration  procedure  was  developed 

using  the  following  principles: 

★  ★ 

-  calculate  and  $2  ^rom  equation  (29)  which  provides  bounds  for 
the  solution  at  the  boundary  (C^$  and  C^$), 

-  calculate  an  estimate  of  $  at  the  boundaries  based  on  4>-j  and 

-  solve  equation  (25)  by  assuming  the  velocity  potentials  at  the 
boundaries.  Since  <t  values  at  the  boundaries  are  specified,  X  values 
are  not  necessary  at  this  step, 

★  ic 

-  using  equations  (38)  and  (39),  calculate  ^  1  and  x2  ^  which  provide 

bounds  for  the  next  estimate  of  x, 

*  *  * 

-  calculate  a  new  estimate  of  x  based  on  x^  ^  and  x^  2, 

-  repeat  the  iteration  procedure. 

The  above  proofs  only  demonstrate  that  at  each  iteration,  a  new  set 
of  upper  and  lower  bounds  are  obtained  for  boundary  fluxes  and  boundary 
velocity  potentials.  The  next  step  is  to  show  that  the  above  iterative 
scheme  is  stable.  Also  it  is  necessary  to  understand  the  important  factors 
effecting  the  rate  of  convergence. 

Computationally,  the  efficiency  of  the  above  scheme  is  based  on  the 
assumption  that  a  large  system  can  be  divided  into  a  series  of  blocks  which 
can  be  individually  stored  in  the  main  storage  area  of  a  computer.  Blocks 
can  be  processed  individually  in  this  scheme  while  the  calculation  of  the 
surface  fluxes  is  done  through  a  relaxation  scheme  for  each  surface  without 
ever  solving  a  large  system  of  equations.  These  computational  considera¬ 
tions  will  later  be  discussed  in  detail. 


D.  The  Convergence  Characteristics  of  the  Iterative  Scheme 

To  understand  the  convergence  characteristics  of  the  iterative  scheme, 
let  us  again  consider  the  simple  two-block  model. 

Start  the  solution  by  assuming  an  initial  flux  vector  at  the  bound¬ 
ary:  A°.  Then,  from  equation  (30),  we  can  calculate  the  potentials 
as  follows: 


*1  =  -l1-!  '  -l1-!,  1-1 


**  =  A"1  p  _  A“V^ 

*2  -2  -1  -2  -2,1-1 


(44) 


•  Calculate  the  potential  vectors  on  the  boundary  surface  and 
average  two  vectors  on  the  boundary  by  using  an  averaging  vector 


a . 


*1,1  =  -1,1  *1 
*2.1  =  -2,1  *2 


(45) 


*B,1  =  a*l,l  +  (1  •  a)  *2,1 

•  Calculate  the  boundary  fluxes  at  the  same  boundary  by  solving  the 
following  pairs  of  equations  for  each  block  by  specifying  poten¬ 
tials  on  the  boundary  surface  or  constraints. 


and 


-1  *l,t  +  -1,1-1 ,1  =  -1 
?1,1  *l,t  =  *B  ,1 


-2  *2,t  +  -2,1  -1,2  =  -2 


(46) 


-2,1  *2,t  =  "  *B,1 


(47) 


In  the  above  equation  values  of  the  right-hand-side  vectors  are  known. 
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★  ★  *  .  J 

•  Average  ^  and  x 2  ^  solutions  at  the  boundaries  with  a  second 
averaging  parameter  8 


x]  =  6X^  +  (1  -  8)  X^  (48) 

By  combining  equations  (44  -  48)  we  can  write  the  following  recur¬ 
rence  relationship: 

Xn+1  =  C  +  DXn  (49) 

where  the  coefficient  matrix  D  can  be  written  as  follows: 

D  =  [bA-j  -  (1  -  8)A2]  [aA'1  -  (1  -  cOA'1]  (50) 

Here  A-|  -j  and  A2  -j  are  reduced  coefficient  matrices  related  to  the  bound¬ 
ary  surface  in  the  following  form: 


-1  -1,1  -1  -1,1 


-2  =  -2,1  -2  -2,1  ( 
The  rate  of  convergence  depends  on  the  eigenvalues  of  the  coefficient 


matrix  D,  which  can  be  simplified  as, 

D=(l  -  a  -  8  +  2a8)  I  -  6(1  -a)^  A2^  -  a(l  -  B )  A^  A^  ^  (52) 

where  I  is  the  identity  matrix.  In  the  case  of  two  equal  size  blocks  with 
A ^  =  A2>  D  can  be  written  as  follows: 

D  =  (1  -  2a  -  28  +  4a8)I 

In  this  case  for  a  =  8  =  0.5,  an  exact  solution  is  obtained  after  one  step. 


E .  Stability  Limits  for  the  Iterative  Solution  Scheme 

To  illustrate  the  bounds  for  convergence  of  the  developed  scheme,  a 
simplified  one-dimensional  model  will  be  utilized.  In  this  case,  the  co¬ 
efficient  matrix  becomes  a  scalar,  d  and  A^A2  matrix  can  be  represented 
by  another  scalar  r.  Of  course  A^A-j,  becomes  1/r.  Then,  the  coefficient 
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matrix  in  equation  (52)  can  be  written  as 

d  =  1  -  a  -  B  +  2aB  -  8(1  -a)r  -  a(l  -  B)/r  (54) 

The  stability  condition  for  this  system,  which  requires 

-  1  <  d  <  1  (55) 

becomes 

2  >  a  +  B  -  2aB  +  B(1  -  a)r  +  a(l  -  B)/r  >  0  .  (56) 

For  example,  a  =  B  =  0.5,  this  condition  can  be  written  as  follows: 

2  >  0.5  +  0.25r  +  0.25/r  >  0  (57) 

Since  r  is  positive  for  positive-definite  systems,  the  stability  condition 
requires 

3-/8<r<3+v^  (58) 


In  a  general  case,  r  represents  the  ratios  of  eigenvalues  of  neighboring 
blocks. 

As  can  be  seen  from  the  above  expression,  the  relaxation  scheme  is 
convergent  only  for  a  narrow  band  of  ratios  of  eigenvalues.  This  means 
that  if  the  neighboring  blocks  are  not  similar  in  terms  of  shape,  size  or 
grid  refinement,  one  cannot  expect  convergence  from  the  iterative  scheme. 

Of  course  in  a  general  system,  a  series  of  eigenvalues  are  involved  which  we 
do  not  have  prior  knowledge  of.  The  above  results  were  obtained  for 
a  =  b  =  O'. 5.  A  more  general  representation  of  the  bands  are  shown  in 
figure  5. 

In  order  to  improve  the  ratio  of  convergence,  we  introduced  a  relaxa¬ 
tion  parameter  w  into  the  scheme  in  the  following  manner: 

•  Once  the  estimates  for  velocity  potentials  were  calculated  and  an 
average  was  computed  in  equation  (45),  the  velocity  potentials  for 
each  surface  vector  was  relaxed  in  the  following  manner: 
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Figure 


(59) 


*1,1  =  “*B,1  +  (1  ‘  *1,1 

*2,1  “  “*8,2  +  (1  -  *2,1 

•  The  second  part  of  equations  (46)  and  (47)  were  then  modified  as 
follows: 

Si  ,1  *l,t  =  *1,1 

-2,1  *2,t  =  "  *2,1  (60) 

•  An  average  for  x  vector  is  again  calculated  from  equation  (48). 

However,  in  this  case,  the  X  vector  is  relaxed  again  with  the  same 

relaxation  factor  to  calculate  the  fluxes  for  each  surface  as  follows: 
1  ^ 

-1,1  =  “-1  +  (1  “  ■)  *1  2  =  “*1  +  ^  1  (61) 

In  this  case  the  coefficient  matrix  D  can  be  written  in  the  following  form: 

D  =  [1  -  ua  -  +  2a>a6]I  -  6w(l  -  a)A^A^^  -  a)a(l  -  gjA^A-j^  (62) 

The  stability  condition  for  the  one-dimensional  case  requires  that 

2  >  a)  [a  +  6  -  2a6  +  6(1  -  a)r  +  a(l  -  6)/r]  >  0  (63) 

In  comparison  with  the  stability  condition  in  equation  (56),  in  this  case 
w  becomes  the  controlling  factor.  Figure  6  illustrates  the  bounds  for 
stability  with  different  relaxation  parameters  w  for  a  fixed  value  of 
a  =  6  =  0.5.  In  the  actual  computations,  as  it  will  be  discussed  later, 
both  were  fixed  at  0.5  and  the  stability  was  controlled  by  only  changing 

iii. 

For  a  general  problem,  r  becomes  the  ratio  of  the  largest  eigenvalue 
for  one  block  to  the  smallest  eigenvalue  of  the  neighboring  blocks  or  vice 
versa.  The  magnitude  of  the  largest  eigenvalue  of  the  coefficient  matrix 
for  a  block  depends  on  the  size  of  the  smallest  element.  On  the  other 
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hand,  the  smallest  eigenvalue  becomes  lower  with  the  increasing  size 
of  a  block.  Thus,  we  expect  small  relocation  parameters  to  be  employed 
when  large  blocks  are  connected  with  small  blocks  with  excessive  grid 
refinements. 

It  remains  an  interesting  problem  to  experiment  with  different  values 
a  and  6  for  each  block  interface  since  there  is  no  requirement  for  a  and 
6  to  retain  constant  values  during  the  iteration  scheme. 
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III.  COMPUTATIONAL  PROCEDURE 
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In  the  above  discussions,  we  tried  to  summarize  the  fundamental 
principles  and  objectives  in  designing  the  present  computational  scheme. 

In  this  section,  we  will  try  to  describe  the  numerical  procedure  in  detail, 
emphasizing  its  computational  aspects.  A  flow-chart  summarizing  the  iter¬ 
ative  scheme  is  shown  in  figure  7.  In  this  flow-chart,  block  surface 
means  one  surface  of  a  block  attached  to  another  single  neighboring  block. 

A  simple  example  for  a  two-block  structure  is  shown  in  Figure  8. 

As  can  be  seen  from  this  flow-chart,  the  first  step  is  to  initialize 
the  density  of  each  element  in  all  blocks  and  assign  a  boundary  flux  dis¬ 
tribution  for  all  global  boundaries  and  inter-block  boundaries  as  defined 
in  equations  (11-12).  In  actual  computations,  we  may  initialize  the  problem, 
for  example,  by  setting  the  densities  to  zero  everywhere  and  by  assuming 
the  flux  boundary  conditions  at  the  upstream  to  be  valid  for  all  block 
boundaries. 

We  then  perform  the  first  set  of  operations.  We  consider  each  of  the 
blocks  one  at  a  time.  There  is  no  order  to  the  sequence  in  which  the 
blocks  are  operated  on.  Given  a  series  of  processors,  these  block  opera¬ 
tions  can  be  distributed  among  these  processors.  For  each  block  we  per¬ 
formed  what  we  called  operation  A,  as  described  in  figure  9.  This  oper¬ 
ation  involved  the  solution  of  the  conservation  of  mass  equation  for  each 
block  individually  under  the  set  of  given  boundary  fluxes  as  specified  in 
equations  (11-12). 

One  problem  in  this  operation,  which  was  mentioned  previously,  was 
the  necessity  to  specify  the  velocity  potential  at  one  grid  point  due  to 
the  singularity  of  the  system.  It  was  pointed  out  that  this  value  could 
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Flowchart  of  the  Block-Structured  Solution  Scheme. 
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be  arbitrarily  chosen  if  the  conservation  of  mass  for  the  entire  block  was 
satisfied  as  stated  in  equation  (14).  For  this  purpose,  we  first  scaled 
the  mass  fluxes  over  the  entire  boundary  which  was  called  operation  C  in 
the  flow-chart  as  shown  in  figure  10.  Details  of  this  scaling  process 
will  be  discussed  later.  After  the  operation  A  was  performed,  a  set  of 
velocity  potentials  was  obtained  and  stored  for  each  block  ($*  in  equation 
(31))  as  individual  sets. 

The  second  set  of  operations  included  the  treatment  of  each  of  the 
surfaces.  This  was  called  operation  E  in  the  flow-chart  and  is  illus¬ 
trated  in  figure  11.  In  this  case,  for  each  surface,  velocity  potentials, 
calculated  from  individual  blocks,  were  picked,  averaged  by  using  equation 
(46)  and  relaxed  as  defined  in  equation  (60).  Again  each  surface  was  ad¬ 
dressed  and  processed  individually.  However,  a  certain  restriction  exist¬ 
ed  here  which  was  due  to  the  fact  that  a  block  may  have  several  surfaces. 

When  processing  a  particular  surface,  data  for  a  particular  block  may  not 
be  readily  available  in  a  parallel  processing  environment  since  it  could 
be  addressed  at  the  same  time  by  another  surface. 

The  last  loop  was  again  a  block  operation.  This  was  called  operation 
B,  as  shown  in  figure  12.  For  each  block,  we  picked  the  averaged  and 
relaxed  values  of  surface  potentials.  We  imposed  these  values  on  bound¬ 
ary  conditions  as  given  in  equations  (39)  and  solved  the  conservation  of 
mass  equation  one  more  time  but  with  this  new  boundary  condition.  The 
outcome  of  this  operation  was  the  boundary  fluxes  for  each  block.  This 
information  was  written  on  a  surface  file  again  to  be  averaged  and  relaxed. 
Thus,  each  time  a  block  was  processed,  its  contribution  to  the  flux  balancinq 
(operation  D,  shown  in  figure  13)  was  simultaneously  performed  as  indicated 
in  the  flow-chart. 
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The  details  of  these  operations  are  listed  below. 


A.  Operation  A 

Operation  A  involves  the  iterative  solution  of  the  potential  flow 
equation  for  each  block  as  defined  in  equations  (10-12).  We  employ  an 
iterative  procedure  for  each  block  as  summarized  below: 

At  each  iteration,  equation  (18)  has  to  be  solved  to  determine  the 
nodal  values  of  the  velocity  potentials.  Rather  than  solving  equation 


(18)  directly,  we  employ  the  following  iterative  procedure: 


.o.n+1  _  r  rt  ,n  «n  n 

Mi  -  Ei  -  -  Mi 


(64) 


Here  n  indicates  the  iteration  step,  where  $  and  x  vectors  are  calculated 
from  the  previous  iteration.  Coefficient  matrix  A^  is  determined  by  using 
the  densities  obtained  from  the  previous  iteration  step.  A?  is  calculated 
at  the  first  step  by  employing  a  fixed  density  distribution. 

The  right-hand-side  vector  is  calculated  for  each  element  of  a  block 
and  assembled  as  a  block  residual  vector.  Matrix  A?  is  never  assembled 
but  calculated  element-by-element  to  be  multiplied  by  a  vector  consisting 
of  the  nodal  values  of  element  velocity  potentials.  The  left-hand-side 
coefficient  matrix  is  a  symmetric,  positive-definite  matrix  which  can  be 
decomposed  only  once  and  stored.  In  this  case,  each  time  a  residual  vector 
is  determined  for  a  vector,  a  forward-backward  substitution  operation  is 
performed  on  this  coefficient  matrix.  The  size  of  the  coefficient  matrix 
is  equal  to  the  number  of  the  grid  points  in  a  block  times  the  half  of  the 
wavefront.  For  example,  for  a  regular  grid  of  10  x  20  x  30  grid  points, 
the  wavefront  would  be  equal  to  10  x  20  x  2  =  80.  If  sufficient  core  is 
not  available  to  store  such  information,  the  coefficient  matrix  has  to  be 
assembled  and  decomposed  at  each  iteration  step. 
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In  our  computations,  considerable  time  was  spent  in  terms  of  formulat¬ 
ing  the  residual  vector.  This  computational  time  can  be  reduced  consider¬ 
ably  if  the  matrices  in  equations  (19-22)  can  be  stored  and  read-in  at 
every  iteration  step.  Only  the  density  has  to  be  recalculated  at  each 
iteration.  In  this  case,  the  core  requirements  increase  by  approximately 
84  x  number  of  elements  for  each  block. 

The  details  of  the  above  iteration  scheme  have  been  previously  presented 
in  detail  [12].  Its  convergence  rate  depends  strongly  on  the  manner  arti¬ 
ficial  density  is  included  in  the  computations.  The  same  scheme  has  been 
applied  without  any  modification  to  the  block-structured  scheme.  As  it 
will  be  shown  later,  no  significant  reduction  in  the  rate  of  convergence 
was  observed  after  the  application  of  the  block-structured  solution  scheme 
as  compared  to  solving  the  same  problem  as  a  single  block. 

B.  Operation  B 

Operation  B  involves  the  solution  of  conservation  of  mass  equation 
for  a  slightly  different  problem.  In  this  case,  the  velocity  potentials 
at  the  inter-block  boundaries  are  specified  and  the  normal  mass  fluxes  at 
the  same  boundaries  have  to  be  determined.  This  problem  can  be  defined  as 
follows:  Consider  the  solution  of  the  conservation  of  mass  equation, 

v*(pV$)  =0  in  n 

with 

pu-n  =  f  on  r° 

and 

=  p  on  rc  (65) 

The  variational  problem  can  then  be  written  as  follows: 

tt  =  dn  +  h°  *ifi  dr  +  /rc  Wp)  dr 


(66) 


In  this  case,  the  Lagrange  multiplier  on  the  boundary  surface  rc  indi¬ 
cates  the  unknown  boundary  flux  to  be  determined  which  corresponds  to  the 
known  velocity  potential  distribution  as  defined  by  p(s,t). 

After  substituting  the  finite  element  approximation,  the  governing 
equations  can  be  written  as  follows: 


Here,  j.  .vector  represents  the  unknown  mass  fluxes  on  each  surface  grid 
i  >  J 

point  that  is  connected  to  a  neighboring  block  surface.  The  solution  of 


equation  (67)  can  be  written  as  follows: 


where  i  is  the  block  number  and  j  is  the  surface  number. 

As  can  be  seen  from  equation  (68),  the  mass  fluxes  at  the  inter-block 
boundaries  can  be  written  in  terms  of  mass  flux  vector  at  global  boundaries 
(F.)  and  the  known  values  of  the  velocity  potential  vector  at  the  inter¬ 
block  boundaries  ( ) -  The  first  part  is  a  constant  vector  through  the 
iterations  while,  d  .  vector  is  recalculated  at  each  iteration.  The  co- 

vl 

efficient  matrix,  in  equation  (68),  is  simply  the  reduced  form  of  the 
matrix  A..  where  all  of  the  grid  points  on  the  boundaries  are  eliminated. 

This  coefficient  matrix  is  different  then  the  one  employed  in  Operation 
A.  Again  by  using  a  constant  coefficient  matrix  scheme,  D  can  be  calculated 
only  once  and  stored.  For  a  block  of  10  x  20  x  30  grid  points  it  will 
have  a  size  of  8  x  18  x  28.  It  will  also  have  a  wavefront  of  2  x  8  x  18. 


There  is  a  basic  difficulty  in  the  above  computational  scheme.  Let 
us  assume  that  a  block  has  six  surfaces  and  all  of  these  six  surfaces  are 
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attached  to  neighboring  blocks.  In  this  case,  the  d.  vector  has  to  be  cal- 

J 

culated  from  all  six  surfaces.  However,  there  is  no  restriction  to  ensure 
that  the  velocity  potentials  calculated  from  each  surface  are  continuous 
along  the  edges  connecting  such  surfaces.  A  unique  value  cannot  be 
assigned  to  $  to  form  the  p.  vector  at  these  points.  Thus,  the  above  set 

”  J 

of  computations  in  equation  (68)  has  to  be  performed  for  each  surface 
individually. 

C.  Operation  C 

Consider  a  single  block  with  a  series  of  boundary  fluxes  on  each  sur¬ 
face.  For  a  rectangular  block  this  would  mean  six  individual  surfaces. 

In  the  present  scheme,  all  of  these  fluxes  are  calculated  individually 
through  computations  involving  the  neighboring  blocks.  We  also  know  that 
we  can  solve  the  conservation  of  mass  equation  only  if  it  is  satisfied 
globally  for  the  entire  block.  To  ensure  such  a  condition,  the  total 
mass  flux  is  calculated  from  equations  (11)  and  (12)  as  follows: 

Ir  f  dr  +  J  c  g  dr  =  e  (69) 

0 

Since  the  boundary  conditions  on  the  global  boundary  is  exact,  the  error 
can  be  distributed  to  the  other  boundary  surfaces  in  the  following  manner: 
JpC  | g 1  dr  =  Q 

g  =  9  -  §  i9l  (?o) 

This  process  ensures  that  the  conservation  of  mass  is  always  satisfied 
before  Operation  A  is  performed. 

D.  Operation  0 

Operation  D  involves  the  balancing  of  fluxes  as  calculated  from 
neighboring  blocks.  As  indicated  in  equation  (64),  at  the  end  of 
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Operation  B,  a  set  of  boundary  fluxes  are  calculated  on  each  surface. 

In  general,  we  can  describe  such  surfaces  with  a  set  of  discrete  points. 
Let  us  assume  that  block  A-|  which  is  joined  to  this  particular  surface 
has  m-j  grid  points  while  block  A^  attached  to  the  same  surface  by  grid 
points.  Let  us  also  assume  that  m^fm^^n  and  all  three  sets  do  not  physi¬ 
cally  match  in  space  with  each  other. 

By  using  a  finite  element  procedure,  we  can  employ  a  shape  function 
as  defined  in  equation  (17)  for  each  surface  in  the  following  form: 


9j(s »t)  =  Nk(s,t)Xi  k(k=l ,n) 

=  Nk(s,t)x1ijik(k=l,m,) 

Vj*5’1'  =  iVs-t>Vj,k(k=,>",2>  <7,> 

Where  X.  .  represents  the  nodal  fluxes  on  the  surface  nodes,  X-  .  .  and 

1  »  K  1  *  J  »  K 

x  .  .  the  nodal  fluxes  on  the  surface  nodes  of  the  neighboring  blocks 

*»*  *3 

i  and  m.  Then,  by  following  the  averaging  scheme  in  equation  (49)  and 
the  relaxation  scheme  in  equation  (61),  it  is  possible  to  evaluate 

2j+1  *  ea,,j +  0  -  pJsu.j 

aft  ■  “af1  +  0  - 


g"+-  - 

am,j 


+  (1 


(72) 


In  the  case  of  blocks  with  matching  grid  points  on  the  neighboring 
surfaces,  the  above  operation  is  conducted  directly  on  the  nodal  values 
of  the  X  vector. 

It  should  be  mentioned  that  the  above  operation  is  in  fact  a  block 
based  operation.  For  each  element  of  the  block,  x  vector  is  calculated 
for  each  surface  and  stored  on  a  surface  based  file  according  to  the 
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equation  (72).  After  each  block  is  processed,  the  outcome  is  the  g.  • 
or  gm  j  vector  on  each  surface.  The  computational  cost  of  the  operation 
is  based  on  the  effort  in  accessing  such  a  surface  file.  For  example, 
when  one  block  is  processed,  the  computed  surface  fluxes  has  to  be  written 
on  each  of  the  surfaces  attached  to  that  block.  The  same  surface  will 
then  be  accessed  by  the  two  blocks  it  is  attached  to.  If  a  surface  infor¬ 
mation  file  is  kept  in  main  memory,  which  includes  both  surface  fluxes  and 
geometry  of  the  elements  on  both  side  of  this  surface,  all  of  the  above 
information  exchange  can  be  directly  made  in  the  main  memory.  Otherwise, 
an  auxiliary  file  has  to  be  accessed  for  storing  such  information.  Know¬ 
ing  the  connectivity  of  the  blocks  and  the  sequence  in  which  this  informa¬ 
tion  has  to  be  accessed,  one  can  optimize  such  an  operation  for  directing 
large  systems. 

E.  Operation  E 

Operation  E  is  very  similar  to  Operation  D  but  somewhat  more  direct. 
Velocity  potentials  are  calculated  from  each  surface,  averaged  and  relaxed. 
Again  a  block  based  information  file  is  generated  which  involves  the 
velocity  potentials  for  each  of  the  grid  points  for  a  block.  At  each  step, 
the  surface  potentials  are  updated  once. 

F.  File  Structure  of  the  Operation 

Two  types  of  basic  information  are  necessary  for  the  operation:  Velo¬ 
city  potentials  and  boundary  fluxes.  This  information  is  stored  for  both 
Operations  C  and  D  in  the  following  manner: 

1 .  Block  file 

This  file  includes  the  following  information  for  each  block: 
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-  Velocity  potentials  for  each  node 

-  Geometry  of  all  the  elements 

-  Shape  functions  for  all  of  the  elements 

-  Element  coefficient  matrices  in  equations  (19  -  22) 

-  Decomposed  Coefficient  matrices  for  a  block. 

This  file  is  read  in  for  one  block  at  a  time.  All  blocks  can  be 

processed  sequentially  or  in  any  order. 

2.  Surface  file 

This  file  includes  the  following  information  for  each  surface: 

-  Surface  fluxes  for  elements  on  both  sides  of  the  surface 

-  Geometry  of  the  elements  on  both  sides  of  the  surface 

-  Shape  functions  for  these  elements 

-  Element  matrices  for  these  elements. 

If  we  review  the  flow-chart,  the  following  I/O  operations  are 
observed:  Operations  A  and  C  require  access  to  the  block  file  while 
Operations  B  and  D  require  access  to  the  surface  file. 

In  the  application  of  the  method  on  an  IBM  4381  conih.  .er  with  12 
megabytes  memory,  the  block  and  surface  files  were  split  into  smaller 
files.  Velocity  potentials  and  boundary  flux  vectors  were  kept  in  main 
memory  at  all  times,  while  the  other  information  was  read  in  when  a  parti¬ 
cular  block  or  surface  was  processed.  This  way  most  of  the  information 
transfer  to  the  neighboring  blocks  in  terms  of  averaging  and  relaxing 
was  performed  in  the  main  memory. 


IV.  DISCUSSION  OF  RESULTS 


To  illustrate  the  applicability  of  the  developed  numerical  procedure, 
a  series  of  test  cases  was  analyzed.  As  discussed  in  the  previous  sections, 
the  convergence  rate  and  bounds  of  convergence  for  this  scheme  depends 
on  the  size  of  the  blocks  and  the  grid  refinement  in  each  block.  One  can 
expect  variations  on  the  size  of  different  blocks  and  the  necessary  grid 
refinement  for  each  block  in  the  case  of  complex  three-dimensional  flows. 

To  test  the  robustness  of  the  scheme,  blocks  were  chosen  in  some  cases  in 
a  way  which  was  considered  to  be  least  favorable  for  convergence  of  the 
scheme.  The  test  cases  include  several  two-dimensional  cases  analyzed 
using  the  developed  problem  as  well  as  a  fully  three-dimensional  case  as 
follows: 

A.  Incompressible  flow  around  a  NACA0012  airfoil 

B.  Compressible,  subsonic  flow  around  a  NACA0012  airfoil  (M^  =  0.6) 

C.  Transonic  flow  around  a  NACA0012  airfoil  (Mro  =  0.83) 

D.  Transonic  flow  around  a  wing-body  configuration  (Mro  =  0.90, 
a  =  3°) 

A.  Analysis  of  Incompressible  Flow  Around  a  NACA0012  Airfoil  (Min  =  0.001) 

Two  types  of  block  structures  were  employed  for  the  analysis  of  this 
test  case.  The  first  case  involved  a  one-dimensional  block  distribution 
in  the  flow  direction  with  three  blocks  and  two  intermediate  surfaces. 

This  block  structure  is  shown  in  Figure  14a.  The  number  of  grid  points 
employed  in  each  of  the  blocks  are  also  shown  in  this  figure.  The  generated 
grid  for  this  block  structure  can  be  seen  in  Figures  14b  -  14d.  The  second 
type  of  block  structure,  shown  in  Figure  15,  involved  twelve  blocks  and 


Figure  14d.  A  section  of  the  3-block  grid  for  the  analysis 
of  flow  around  a  NACA0012  profile. 


No.  of  Grid  Points 


Figure  15.  Twelve-block  grid  structure  for  the  analysis 
of  flow  around  a  NACA0012  profile. 


twelve  intermediate  surfaces.  In  this  case,  blocks  of  different  sizes 
were  placed  in  all  three  principal  directions  to  test  the  robustness  of 
the  numerical  scheme.  The  assembled  grids  for  both  types  of  block  struc¬ 
tures  are  identical  in  both  cases  and  are  shown  in  Figure  16.  Both  grids 
were  considerably  coarse  with  only  ten  elements  over  the  airfoil  but  they 
still  provided  adequate  accuracy  for  the  solution  of  this  simple  program. 

The  numerical  results  obtained  for  an  inlet  profile  of  M.  =  0.001 

in 

from  the  three-block  and  twelve-block  grid  structures  were  compared  with 
the  exact  solution  of  this  problem,  as  shown  in  Figure  17.  In  both  cases 
no  relaxation  was  employed  and  both  required  about  35  iteration  steps  for 
convergence. 

B.  Analysis  of  Compressible  Subsonic  Flow  Around  a  NACA0012  Airfoil  (Min  =0.5) 

For  the  solution  of  this  test  problem  the  same  three-block  and  twelve- 
block  grid  structures  described  in  case  A  were  employed.  The  obtained 
numerical  results  are  compared  in  Figure  18  with  results  from  a  single¬ 
block  C-type  grid  with  (61)  x  (50)  grid  points. 

The  agreement  with  the  one-block  solution  is  very  good  considering 
the  coarseness  of  the  grid.  The  number  of  iterations  required  for  conver¬ 
gence  was  about  40  for  both  cases  and  no  relaxation  factor  was  employed. 

C.  Analysis  of  Transonic  Flow  Around  a  NACA0012  Airfoil  (M^  =  0.83) 

Again  two  types  of  block  structures  were  employed  in  the  analysis. 

First,  a  three-block  structure  shown  in  Figure  14a was  tested  with  blocks 
distributed  only  in  the  flow  direction.  The  second  case,  shown  in  Figure 
19,  was  a  twelve-block  structure  involving  blocks  distributed  in  two  dimen¬ 
sions.  This  time  two  blocks  were  employed  over  the  airfoil  surface  to 
provide  a  complicated  test  case  for  locating  the  shock.  The  grids  pro- 
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Figure  19.  Twelve-block  grid  structure  for  the  analysis  of 

transonic  flow  (M^n  =  0.83)  around  a  NACA0012  profile. 
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duced  for  both  block  structures  were  identical  but  considerably  finer  than 
those  employed  in  the  ones  used  for  incompressible  flow  problems.  A 
cross-section  of  the  grid  is  shown  in  Figure  20a.  In  Figure  20b  the  grid 
distribution  over  the  airfoil  surface  can  be  seen.  The  grid  employed 
( 71 )x( 16x( 3)  =  3308  total  grid  points  with  a  30  element  distribution  over 
the  airfoil.  The  distribution  of  elements  for  each  block  is  shown  in 
Figures  14a  and  19. 

The  numerical  results  obtained  from  the  three-block  and  twelve-block 
grid  structures  were  compared  with  results  obtained  with  a  single-block 
C-type  grid  and  are  shown  in  Figure  21.  Both  results  compare  favorably 
with  the  single-block  solution,  obtained  from  the  C-type  grid  shown  in 
Figure  22.  The  number  of  iterations  required  for  convergence  was  about 
450  in  both  cases  and  involved  a  variable  relaxation  factor  to  in  the  range 
of  .01-. 55. 

D.  Transonic  Flow  Around  a  Wing-body  Configuration  (M^  =  0.90,  a  =  3°) 

A  24-block  grid  structure  was  chosen  for  analyzing  the  transonic  flow 
problem  around  a  wing-body  configuration  as  shown  in  Figure  23.  The 
details  of  the  generated  grid  are  shown  in  Figure  24.  The  number  of  grid 
points  in  each  block  is  listed  below: 


Block  No. 

No.  of  Grid  Points 

Bill 

7x8x4 

B121 

14x8x4 

B131 

7x8x4 

B112 

7x8x4 

B122 

14x8x4 

B132 

7x8x4 

B211 

7x8x1 

B221 

14x8x1 

B231 

7x8x1 

B212 

7x8x1 

B222 

14x8x1 

B232 

7x8x1 
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Figure  23a.  Block  structure  for  analyzing  the  transonic 
flow  around  the  wing-body  configuration. 
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stribution  for  a  single  block  with  and 
hidden  lines  (Block  #B422  from  Figure  23b). 


Expanded  view  of  the  block  grids  for  the  first 
column  of  blocks  for  the  wing-body  problem. 
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Figure  24g.  Expanded  view  of  the  block  grids  for  the  third 
column  of  blocks  for  the  wing-body  problem. 


Block  No. 

No.  of  Grid  Points 

B311 

7x8x6 

B321 

14x8x6 

B331 

7x8x6 

B312 

7x8x6 

B322 

14x8x6 

B332 

7x8x6 

B411 

7x8x4 

B421 

14x8x4 

B431 

7x8x4 

B412 

7x8x4 

B422 

14x8x4 

B432 

7x8x4 

The  total  number  of  grid  points  is  8207.  Comparison  of  obtained  results 
as  obtained  from  the  previous  solution  of  the  same  grid  as  a  single  block 
[6]  is  shown  in  Figure  25. 

As  can  be  seen  from  these  results,  the  24-block  solution  reproduces 
the  single  block  solution  reasonably  well.  The  differences  between  the 
two  results  can  be  attributed  to  the  coarseness  of  the  grid  in  terms  of 
computing  mass  fluxes  between  the  blocks. 

The  number  of  iterations  required  for  convergence  in  this  case  was 
550  and  involved  a  variable  relaxation  factor  u  in  the  range  of  .01 -.75. 


Comparison  of  the  pressure  distribution  at  the 
wing  root  fo,r  single  block,  multi-block  and 
experimental  results. 


V.  CONCLUSIONS 


The  procedure  described  in  this  report  involves  the  solution  of  tran¬ 
sonic  potential  flow  equations  using  a  block-structured  approach.  The 
basic  objective  of  this  research  was  to  develop  methodology  for  analyzing 
large  flow  problems  involving  complex  geometries.  The  developed  capabili¬ 
ties  include, 

•  the  generation  of  a  computational  grid  around  a  complex  geometry 
based  on  a  block-structured  grid  generation  scheme, 

•  the  analysis  of  such  a  system  by  using  a  block-structured  solution 
scheme  for  transonic  potential  flows. 

The  grid  generation  procedure  was  presented  in  previous  publications. 

The  block-structured  solution  scheme  is  presented  in  this  report  only  for 
potential  flows.  It  involves  the  solution  of  conservation  of  mass  equation 
on  a  block-structured  basis  in  terms  of  the  velocity  potential  <j>.  Each 
block  is  stored  independently  and  analyzed  in  the  computer.  Iterations 
involve  solution  of  conservation  of  mass  equation  for  each  block  as  well  as 
balancing  the  mass  flux  between  the  neighboring  blocks. 

The  procedure  is  expandable  to  the  solutions  of  Euler  and  Navier- 
Stokes  equations.  In  the  case  of  Euler  equations,  equation  for  the  conser¬ 
vation  of  entropy  will  also  have  to  be  solved.  This  again  will  involve  the 
solution  of  the  equation  for  each  block  as  well  as  balancing  the  entropy 
to  the  neighboring  blocks  in  the  flow  direction.  Such  an  approach  provides 
the  capability  of  solving  different  types  of  equations  approximating  the 
flow  field  on  the  same  grid  but  for  different  blocks.  For  a  complex  air¬ 
craft  configuration,  the  capability  to  solve  potential,  Euler  or  Navier- 
Stokes  equations  at  different  flow  regions  is  necessary  considering  the 
present  computer  hardware  capabilities  and  also  the  complexity  of  the  problem. 


As  stated  above,  the  potential  flow  solution  of  transonic  flows 
requires  the  only  the  treatment  of  the  conservation  of  mass  equation.  In 
terms  of  solving  Euler  or  Navier-Stokes  equations,  the  same  equation  has  to 
be  treated.  In  fact,  the  elliptic  nature  of  the  conservation  of  mass 
equation  governs  the  convergence  rate  of  relaxation  schemes  for  all  three 
types  of  formulations.  Therefore,  the  efficiency  of  the  computational 
scheme  in  terms  of  solving  the  conservation  of  mass  equation  for  large 
systems  is  of  general  importance  for  solving  Euler  and  Navier-Stokes 
equations. 

Currently  the  method  is  being  extended  to  the  solution  of  Euler 
equations  including  solving  coupled  potential  and  Euler  solutions  at 
different  flow  regions.  Another  important  consideration  is  the  application 
of  the  numerical  procedure  in  a  parallel  processing  environment.  As 
described  in  this  report,  the  processing  of  blocks  and  surfaces  can  be 
performed  in  any  order.  This  provides  considerable  flexibility  in  defin¬ 
ing  an  optimum  solution  procedure  for  a  given  parallel  processor  hardware 
configuration. 
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